home *** CD-ROM | disk | FTP | other *** search
-
- // This function solves Kepler's equation by the Newton
- // method.
- // Inputs: m -- mean anomaly
- // e0 -- starting estimate for eccen anom
- // e -- orbit eccentricity
- // mi -- maximum number of iterations
- // t -- stopping tolerance
- // Outputs: sk returns the final estimate for eccen anom
- // * Input m and e0 and output sk are in degrees. *
-
- double sk(double m,double e0,double e,int mi,double t)
- {
- /* Solve Kepler's equation by the Newton method. */
- double en,ep;
- int i;
-
- e0 = e0 * DTR;
- m = m * DTR;
-
- i = 0;
- en = e0;
- do {
- i++;
- ep = en;
-
- // The next line of code is the Newton iteration:
- // x(n+1) = x(n) - (1/(dy/dx)) * y(n).
-
- en = en - (en - e*sin(en) - m)/(1 - e*cos(en)) ;
- en = fmod(en,TWOPI);
- if ( fabs(en-ep) < t) break;
- } while (i <= mi );
- return(en*RTD);
- }
-
-
- SAMPLE INPUT/OUTPUT
- Initial eccentric anom = 0
- Input mean anom = 286.879298
- Input maximum its = 15
- Input stop tolerance = 0.000001
- Input eccentricity = 0.100000
- Final eccentric anom = 281.260008
- Mean anom computed from solution = 286.879298
-
- This solution took five iterations; the solution E was
- checked by plugging it back into Kepler's equation
- and computing M.
-